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We study the distribution of vortex Hues in a three-dimensional lattice with inhomogeneous cou- 
plings. We investigate the spatial distribution of the number of vortex lines, showing how the 
vortex lines are expelled from the region with higher couplings. Results for the correlation functions 
are presented, together with a discussion of the implications of the obtained results for the vortex 
distribution in ultracold trapped gases and in neutron stars. 

INTRODUCTION 

The study of vortex and vortex lines in superfluids and superconductors is an evergreen field of research [iHSI . In 
type-I superconductors, the rigidity of the vifavefunction of superconducting electrons determines the expulsion of the 
magnetic field (Meissner effect): the superconductivity is destroyed at a critical magnetic field. However, when the 
value of the Ginzburg-Landau parameter exceeds a critical value, then the magnetic field can enter the superconductor: 
the flux penetrates in a regular array of flux tubes, each carrying a flux quantum HUHj. The corresponding properties of 
vortices in type-II superconductors have been then intensively studied in BCS and high-temperature superconductors 

[SHS] 

Vortices have been extensively studied also in neutral superfluids: as it is well known, the effect of a rotation on 
neutral particles is equivalent to that of a magnetic field on a system of charged particles [3]. As a consequence, one has 
in neutral superfluids the counterpart of the effects regarding the penetration of a magnetic field in a superconductor: 
the expulsion of the magnetic field has its equivalent in that a neutral superfluid, up to a certain value of the rotation, 
does not rotate together with a rotating trap or vessel. Thus the use of Ginzburg-Landau theory to describe vortex 
states is not restricted to superconductivity and it can be applied, with certain adjustments, to neutral superfluids. 

The study of vortex states in neutral superfluids witnessed in the last two decades a substantial additional interest 
due to progress in the realization, control and manipulation of ultracold bosonic and fermionic gases [HI . Corre- 
spondingly, the study of properties of vortices in ultracold gases subjected to rotating external potentials attracted 
a significant attention [9l-fTT]. For trapped Bose-Einstein condensates, when the condensate rotates with angular 
velocity Q higher than a critical value, one or several vortices nucleate: for more rapid rotations, the vortices form a 
dense triangular array [11) . At even higher rotating frequencies the superfiuid properties are lost and when Q is close 
to the radial trap frequency the lowest-Landau-level approximation becomes applicable [H]: the rotating dilute gas 
is expected to exhibit a transition from a superfiuid to various strongly correlated, non-superfluid states analogous to 
those typical in the fractional quantum Hall effect for electrons in a strong perpendicular magnetic field |11H13) . 

We also mention that the recent experimental realization of synthetic magnetic and electric fields acting on neutral 
atoms using a suitable spatial variation and time dependence of the effective vector potential 1141 115] opens the way to 
simulate many-body systems in controllable (static) electro-magnetic fields [TB] : such a possibility, eventually together 
with the use of rotating traps [HI [T^, opens the way to the study of bosonic and/or fermionic systems in a vast class 
of simulated magnetic fields |17| . including spin-orbit coupling. 

The study of vortices has been carried out both in two- and in three- dimensional systems. In two dimensions, the 
thermally generated vortices form, at low temperature, pairs of opposite vorticity, giving rise to power-law correlation 
functions. At the critical point of the corresponding Berezinskii-Kosterlitz-Thouless (BKT) transition the vortices 
start to unbind and high-temperature exponential decaying correlations arise |18l I19| . A paradigmatic model in which 
the BKT transition has been deeply investigated is the XY model on a lattice, which is also commonly used to describe 
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defect-mediated phase transitions and superfluidity in two dimensions [20] . The unbinding of the vortices around the 
XY transition point was studied and visualized using Monte Carlo simulations [21]. The XY model description has 
the advantage that it is possible to extract the topological degrees of freedom (i.e., the vortices) by constructing the 
dual theory [52] , showing that the vortices interact via a logarithmic potential. A physical realization of the lattice 
XY model is provided by Josephson junction arrays, i.e. superconducting grains located at the vertices of a network 
[531 [23]. Furthermore, for large number of particles per site and negligible fluctuations of the number, the effective 
Hamiltonian of ultracold bosons in optical lattices can be mapped to that of an (eventually quantum) XY model 
[25H27j : the vortex number and the BKT transition were experimentally detected in a two-dimensional Bose gas in 
an optical lattice at finite temperature [55] . 

Passing from two to three dimensions, e.g. by coupling two-dimensional systems, then vortex lines connect the 
core of the two-dimensional vortices: the vortex lines must either close on themselves (vortex rings) or terminate on 
boundaries [S]. Following the original ideas put forward by Onsager [53] and Feynman [30] that vortex excitations 
are responsible for the superfluid transition in '^He, the properties of vortices in three-dimensional superfluid systems 
were subsequently investigated: a real-space renormalization technique was used in |3lj to generate screened vortex 
energies and core sizes. The XY model is useful also to study vortex excitations in three dimensions: using the 
dual description, it is possible to write the Hamiltonian in terms of vortex loop degrees of freedom [55]. Monte 
Carlo results highlighted the role of vortex loops, showing indeed that they are responsible for the phase transition 
in the three-dimensional XY model [32]. Using a three-dimensional generalization of the two-dimensional Kosterlitz 
real-space scaling procedure, Shenoy wrote down scaling equations for the vortex-loop fugacity and the loop coupling 
[33]. This computation was then generalized to the anisotropic three-dimensional XY model |34j in order to study 
the 3D-2D crossover. 

Important contributions to theory of vortices and vortex lines were made by Luciano Reatto and his collaborators 
[35H38j : in [35j the energy of a ^_ffe vortex line and the density distribution in the core region were computed using 
model many-body wavef unctions. The velocity of translation of the vortex ring was also computed and shown to 
be in agreement with experimental results. The variational energy was calculated using the Percus-Yevick and the 
convoluted-hypernetted-chain integral equations: the structure factor entering these equations was taken by measured 
values and assumed linear, with slope hk/2mc, at small wavevectors |35| — as it is well known, a model wavefunction 
leading to this result was constructed in the seminal paper "39]. The properties of vortex lines in *He were reconsidered 
in |371 138| where a shadow wavefunction was used (optimized shadow wavefunctions for the ground state of liquid and 
solid ^7?e and at the liquid-solid *He interface were studied in [40H42] ). In [38] the excitation energy of the vortex 
was then computed showing a very weak dependence on the density when the latter is increased from equilibrium to 
near freezing, giving insights on the properties of the classical core parameter [38 . 

The problem we address in the present paper is the study of the effects of inhomogeneities on the vortex lines 
distribution in three-dimensional lattices. As motivation of our work, we have in mind mostly two systems: trapped 
ultracold atoms in optical lattices and neutron stars. Regarding the former, the low-energy properties of ultracold 
atoms in deep optical lattices are in general described by Hubbard-like models both for bosons [33] and fermions 
[33]- When the system is superfluid, in each well i of the periodic potential one can define a phase 0^: either due to 
temperature and/or to the effect of a rotation, vortex lines will develop in the system. Since an external potential 
is usually present in such trapped systems [3 llOj , then an arising issue is the determination of how the trapping 
potential affects the distribution vortex lines. 

The study of the effect of inhomogeneities on the vortex lines in superfluids in a lattice is also relevant in the context 
of neutron stars (NS). NS are compact objects, with typical masses of order 1.5 solar masses and typical radii of order 
10 km, formed from the collapsed core of fairly high mass stars when they explode as supernovae. Apart from their 
extreme compactness, observed NS are usually characterized by very high rotation frequencies — up to a kHz — and 
very high magnetic fields — up to lO^'^-lO^^ G. A NS can be roughly subdivided into three zones having different 
composition and properties: the first is a shell called the outer crust composed of a lattice of atomic nuclei permeated 
by an electron gas. The outer crust ends at about 100 m from the surface, where the density reaches 10^^ gcmT^ . The 
next layer in is called the inner crust and consists of a lattice of very neutron rich nuclei permeated by free neutrons 
and electrons. Theoretical calculations [33 |3S] and indirect observational evidence [33 |3H] support the idea that 
crustal free neutrons form a superfluid. The space dependence of the superfluid gap has been also studied I46| . 
The crust ends at about 1 km from the surface where the density reaches 2.7 x 10^"* gcm^'^ = 0.16 /m~'^. The inner 
region, called the core, is formed by free neutrons, protons and leptons (neglecting any more exotic composition). 
Again observations seem to suggest that neutrons are superfluid while protons are superconducting [471 148] . 

Vortices are expected to appear in both the inner crust and the core of a NS. Moreover, they are expected to 
interact in the crust with normal phase lattice nuclei [491 . This microscopical interaction is widely believed |50j to be 
the mechanism that gives rise to glitches. A glitch is a sudden increase of the otherwise slowly and steady declining 
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neutron star spin period (for a recent statistical study of glitches see [5T]). The inhomogeneous neutron superfluid fills 
a layer spanning three order of magnitude in density and it feels an effective periodic potential created by the crustal 
lattice: free neutrons interact with the lattice nuclei feeling a one-body potential that can be thought as having a 
Thomas-Fermi like shape [52] . At different densities both the lattice nuclei and the lattice length are different [S3J [S3] : 
in any case the potential becomes negligible (because of the fine range of the nuclear force) without superimposing 
with those of the first neighbours, creating in this way a (deep) lattice potential for the superfluid neutrons. An 
interesting issue is then the study of the effect of the inhomogeneity on the vortex lines in the inner crust of a rotating 
NS and possibly the individuation of a simplified model capturing the main properties of vortex lines in NS. As 
a first step towards the individuation of such a simplified model (having parameters derived from the microscopic 
Hamiltonian), we study in the following the distribution of thermally excited vortex lines for a superfluid in a lattice 
with inhomogeneous couplings proportional to the superfluid gap in the NS inner crust: for this purpose we use the 
data for the gap presented in Ref. |i5] . 

To explore the main properties of vortex lines in inhomogeneous lattices and have results independent of the micro- 
scopic details, we decided to investigate the inhomogeneous version of the XY model described by the Hamiltonian 

H = - ^Jij cos (f>j). (1) 

In Eq.Q, the sum is on the nearest-neighbours pairs of sites z, j of a three-dimensional lattice (supposed cubic for 
simplicity) and is a phase variables, having values between and 27r. Furthermore, Jij is the coupling between the 
sites i and that we will assume to be inhomogeneous in the space (typically gaussian in one or all directions). We 
will mostly consider in the following open boundary conditions. 

Several reasons lead us to consider Hamiltonian ([l]) to model the effect of inhomogeneities on vortex lines: first, 
in the homogeneous limit {Jij = J across the lattice), the XY model can be expressed in terms of the topological 
degrees of freedom in two and three dimensions, giving important insight on vortex-vortex interactions and on the 
behaviour of correlation functions. Moreover, the Hamiltonian ([l]) could be realized on its own by using Josephson 
junction arrays with Josephson energies Jij varying across the lattice or by considering ultracold atoms in an optical 
lattice, whose Hamiltonian can be mapped in a suitable range of parameters to ([I]). As a last reason, more important 
for our purposes, we observe that in the presence of a superfluid in a lattice, the macroscopic wavefunction ipj can 
be written as ipj = y^nje"'^^ : then the Hamiltonian will have a term cx ~2t^^,-^-^ y/^i^j cos {(j>i — (f>j), where t is the 
hopping parameter. This means that we can proceed with the identification 

= 2t^n,nj. (2) 

Then, if we have a superfluid in which the number rii is fixed, e.g. by an external potential, and we want to focus on 
the equilibrium and dynamical properties due to the behaviour the phases, we can stick to the Hamiltonian ([T]) with 
the identiflcation 

The plan of the paper is then the following: in Section | we remind the main properties of the Hamiltonian ([T]) 
in the homogeneous limit ( Jij — J) and we recall how the vortex degrees of freedom are introduced. The effect 
of inhomogeneities is addressed in Section [| where we present numerical results for the spatial correlations and the 
number of vortex lines, pointing out that the vortex lines are expelled from the region with large couplings. Results 
with an external magnetic field are also presented. Our conclusions, together with a discussion of perspectives and 
future work, are presented in Section! 



THE HOMOGENEOUS LIMIT 



In this Section we recall the main properties of the Hamiltonian ([T]) in the homogeneous limit, i.e. when the 
couplings Jij are taken equal to some value J, so that ([T]) reads 

= -J^ cos (0, -</.,). (3) 

We observe (jS) can be written a.s H = -~JJ2{ij) ' ^j^ where Si — (cos 0^, sin 0^). A detailed discussion of several 
properties of the three-dimensional XY model is in |55j . 

The XY Hamiltonian (|3| displays a phase transitions, belonging to the 3-D — XY universality class at a critical 
temperature Tc in which the order parameter (cos0i) becomes non- vanishing: accurate Monte Carlo estimates give 
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[SSI [57] 

keTc ~ 2.202J, (4) 

i.e. PcJ - 0.454 (where Pc ^ l/^s?;)- 

A simple estimate of the T^. can be obtained by mean- field [SS]: the mean-field Hamiltonian reads 

Hmf = — zm J cos (j>i, (5) 

i 

where z = 6 is the number of nearest-neighbours and m = {cos is the parameter to be self-consistently determined. 
One has 



cos e 



fizjrn cos cj} 



it follows 



Jm cos 

Ii {PzJjn) 



Iq {(3zJm) ' 



(6) 



(7) 



where /„ (a:) is the n-th modified Bessel function [59j . By performing the limit m — >■ 0, one finds the mean-field critical 
temperature 

which should be compared with Q. 

A useful way to go beyond mean-field results, both for the critical temperatures and the critical exponents, is 
provided by rewriting Hamiltonian ([s]) in terms of the correct topological excitations, the vortex loops. To this aim, 
one has to resort to dual transformations [22] : details are found in [33] [331 [SOI [H] ; here we sketch the main steps in 
order to show how vortex loops emerge in the effective description of the system. 

The dual transform procedure is the same as in the two-dimensional XY model: the main point consists in 
expanding e~^^ in terms of the basis eigenf unctions invariant with respect to the symmetry of the Hamiltonian |22) . 
ft is convenient to write the Hamiltonian ([s]) as 

iJ = -J^cosA^^j, (9) 

where is the discrete derivative in the direction /i, with /i — x,y,z denoting the versors in the three directions 
(e.g., Ax(f>j — cos(/)j+£ — cos0j, and so on). To write the Hamiltonian ^ in term of the dual variables, which turn 
out to be the vortex lines themselves, one has to introduce the lattice bond variables Uj^^ (which are integer variables) 
and then integrate over the initial variables (pi [55]. Expanding the Boltzmann factor e~^^ in terms of e'^j ni,^^^'t>i 
and integrating over the 0i one gets a constraint, as in any dual transformation [22j : the constraint to be satisfied is 
^„ A^rij^^ = A • n = 0. One has then to write Uj^^ as the curl of a dual lattice bond variable in order to satisfy the 
constraint: using the Poisson summation formula [22] amounts then to introduce the vortex loop variables J [r). This 
procedure is the equivalent of the duality transformation that in 2D defines the Coulomb gas Hamiltonian, where the 
charges of the Coulomb gas are the (2D) vortices [22 • One finally gets for the partition function 

Z = Y,e-0'^: (10) 
{J} 

the effective Hamiltonian is written in terms of the vortex loop variables J (f) and it reads 

where U (r) is the three-dimensional lattice Green's function U (f) (apart from a constant): U (r) = U {r} — U{Q), 



where for large distances it is J7 (r) ^ 1/ | r | [55J. In Eq.(ll I Kq ~ J at low temperatures (the full dependence upon 
temperature is discussed in [55(). 
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The vortex loop variables J(r) are vectors, having the three components (r) (ji = x,y,z). The components 
(f) in each site r of the dual lattice are integer variables 

J^(rO = 0,±l,±2,---: 

at low temperatures only the values (r) — 0,±1 are relevant. One has the following local conservation law: 

A • J (f) = 0, (12) 

meaning nothing but that the loops are closed. The physical meaning of the variables (r) is that a vortex passes 
from one elementary cell cube of center r to a neighbour cube of center r + /i if the sum of the phase differences 
{mod2Tr) across the square face separating the two cells has modulus 27r, i.e. if a vortex line passes through the two 
centers (from r to r + /i if (r) = 1, and from r + /i to f if J^^ (r) = — 1). In the 2D case the vector J (r) is a scalar, 
meaning that one has a counter-clockwise (clockwise) vortex around the plaquette center r if the value is 1 (—1), i.e. 
one finds back the standard definition of 2D vortices. 

An important consequence of rewriting the Hamiltonian in terms of the correct degrees of freedom (r) is that it 



is possible, starting from (111, to perform a real-space renormalization study: scaling equations are obtained for the 



vortex-loop fugacity y{£) and the coupling IC{£) as a function of the scale £ [33]. The final result is 

^=IC-AyJC^; % = (Q-^JCC)y (13) 



where A = 47r^/3 and /! « 1 -I- In [a/ac {£)], with a = and Oc {€) the finite core size at scale £. Analysing Eqs.(13|, 
one can have a good estimate of the critical temperature and of the critical exponents |33j . 

The emerging physical picture is the following: the three-dimensional XY model has vortex loops as topological 
excitations, with loop segments interacting via a 1/r potential. Thermally created loops are few and small at low 
temperatures: when the temperature is increased the loops expand and smaller loops wedge between larger loops. 
This screening weakens the binding of larger loops, until the dominating largest loops are expelled at the critical 
temperature, with a proliferation of vortex loops. 



EFFECT OF THE INHOMOGENEITY & NUMERICAL RESULTS 

In this Section we study the Hamiltonian ([ij) with inhomogeneous couplings Jij: to fix the notation, we consider in 
the following three kinds of inhomogeneity. First, we consider a gaussian radially symmetric coupling given by 

J,^. = Joe-'■^/2^^ (14) 

where r is the distance (hereafter measured in units of the lattice length) from a chosen origin of the midpoint of the 
sites j, j; Jo is an overall coupling and i? is a length quantifying the spatial extent of the couplings Jij. The choice 



( 14 1 refers to the physical situations in which an external potential forces the superfiuid to stay close to the centre, 
as it is the case, e.g., for trapped superfluids in parabolic potentials. 

Another possible choice amounts to considering the inhomogeneity only in one direction: we set 

J,, = Joe-^'/2fl^, (15) 



where z is the distance from a chosen origin along the z direction of the midpoint of i and j. The couplings (15) 
correspond, in an ultracold atom setup, to a trap squeezed in the z direction and refer to a superfiuid occupying a 
certain "slice" of a three-dimensional system: this is the case also for the neutron superfiuid in the crust of a NS. 
For this reason, we also consider as a third choice couplings Jij proportional to the superfiuid gap in the NS inner 
crust: we use the data for the gap in Ref. [45] obtained using the /g correlated model (see table 1 in [45]). Since in 
a NS the critical temperature is order of MeV and the temperature T is much smaller than T^, the relevant regime 
of temperatures to model NS through the Hamiltonian ([T]) is the low-temperature one. Denoting with z the radial 
direction of the NS, the data from [IS] as a function of the distance z from the centre of the inner crust are well fitted 
with a function {A + Bz) e~^'' . We will use then (with Jq = -R-S and zq = RA/Jq) 



J^,=Jo^e-^^-^'>y'/'^\ (16) 
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It should be noted that even if the XY model is usually formulated with two energy scales J and T, only their ratio 
has a physical meaning. Therefore all of the results only depend on the adimensional parameter /3Jo (or jSJ in the 
homogeneous case). 

As a starting point of the following discussion, we observe that in the homogeneous limit [i.e., Jij = J in Eq.([l])] the 
critical temperature Tc is proportional to J: for temperatures smaller than Tc a few small vortex loops are thermally 
generated. If there is a shallow dependence of the couplings, one can consider that there is a space dependent value 
of the critical temperature for the vortex loop blow out. Then, with the couplings varying across the lattice, as a 
consequence of the inhomogeneity, in the regions with small (large) coupling constants there are many (few) large 
(small) vortex loops. This can be rephrased in terms of a local density approximation (LDA), where the effective 
local temperature is T'effl'^) ~ JqT/ J{r)- According to this picture, there should be many vortex loops in the regions 
in which the coupling is very small, e.g., for large r with the couplings ( 14 ): however, there are no superfluid particles 
there, as shown by Eq.(|2|, and an effective number of vortices will reach a maximum value at intermediate values of 
the couplings. One can think to the presence of regions with a maximum effective number of vortices as related to the 
occurrence of vortex-vortex interactions that are no longer translationally invariant and effectively attractive in the 
regions with intermediate values of the couplings. One can quantitatively test the validity of the LDA by computing 
the local observables of interest in the homogeneous case for various temperatures and comparing the results with the 
local values of the same observables in the inhomogeneous system, with coupling /? J(r) as we show later in Fig. [9] 

To quantitatively test and visualize the described scenario involving the effective exclusion of vortex loops from the 
regions with higher couplings, we numerically studied the Hamiltonian ([I]) with the couplings (14)-(16) by means of 
a Monte Carlo simulation with open boundary conditions. Vortex lines were identified using the standard definition 
[55] : a vortex passes from one elementary cell cube to a neighbour cell if the sum of the phase differences {mod 2tt) 
across the square face separating the two cells has modulus 27r. In other words, if the face between two cubes has 
the four sites 1, . . . , 4 and phases 0i, . . . , 04, then one computes the quantity (l/27r) J^a ^a4> where a = 1, . . . , 4 and 
Aq0 = 4>a+i — (mod 27r) (with 05 = (pi). The winding direction of the vortex line depends upon the sign of this 
sum (not shown in the figures). 

We examined the typical vortex configurations by taking snapshots. Figs. [T][7] show probable configurations in 
various conditions using 21^ sites. In Figs. [l][7] the sites of the cubic lattice are depicted as points and three- 
dimensional images of the lattice with the vortex lines are presented: the vortex lines, connecting the centres of the 
elementary cell faces, are drawn as solid lines. In presence of inhomogeneous couplings, the thickness of the lines is 
faded proportionally to the couplings Jij, i.e. the full solid lines correspond to the maximum value of the couplings 
and no line is drawn in absence of coupling: different intensities of grey correspond to intermediate couplings |62j . 

The simulation with homogeneous coupling is shown for comparison and it illustrates the nature of three-dimensional 
XY phases transition occurring at 13c J ^ 0.454, as shown in Figjl] (below T^) and Fig[2] (above Tc). The transition 
point is also signalled by the behaviour of the order parameter m — {cos(f>) (not reported). Long vortex lines such 
as the ones in Fig. [3] are unstable below T^. this configuration has been prepared with a sudden quench in Monte 
Carlo dynamics from very high to very low temperature (/3J quenched from to 10). Notice that such configurations 
are metastahle in the Monte Carlo dynamics and can significantly increase the thermalization time: open lines have 
a very long time scale since they can only annihilate touching the boundaries, even if they are not very sensitive to 
the presence of the boundaries themselves. In the same way as they can survive for many simulation steps, they may 
take a similar amount of time to nucleate, so special care has been taken by varying the temperature slowly in the 
simulations. E.g., when crossing the critical temperature, open lines can form by joining several closed loops, so that 
if the temperature starts slightly below the density of closed loops is high enough to allow rapid thermalization. 
In the same way, a rapid quench from above the critical temperature to a very low temperature does not allow the 
breaking of long vortex lines that can remain trapped in the system for a long time (as shown in Fig. [3|. 

In presence of an inhomogeneous coupling, regions with a higher coupling tend to expel vortices. In Fig. [4] we show 
the effect of the spherically symmetric inhomogeneity ( 14 ) : the central core remains vortex free because the coupling 
is higher. Outside the core the vortex density grows rapidly and reaches a constant. However, as the coupling grows, 
the energy carried by vortices decreases. Eventually the vortex loses its meaning when the coupling goes to zero, as 
the phases 4> are decoupled and only a finite shell gives a contribution to the thermodynamic properties. To better 
illustrate this point, we decided to depict the vortex lines as solid where the coupling is maximum and make them 
fade proportionally to the coupling. 

We also examined the case of a planar geometry defined by the Eq.(15|. The region forbidden to vortices is a slab 
region of size R forming a pancake geometry: we observe the creation of two vortex walls, as illustrated in Fig. [5] 
The inner crust acts as a potential barrier for vortices: if the temperature in the centre is lowered, small vortex loops 
can tunnel through. 

The presence of a strong enough magnetic field can make vortices inside the barrier as shown in Fig. [6j The 
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FIG. 1: A Monte Carlo snapshot in the homogeneous limit for a temperature below the critical temperature (/3J = 0.6, while 
j3cJ — 0.454). The vortices are small in radius and clearly separated. In this figure and in the foUowings the sites of the cubic 
lattice are depicted as points in a three-dimensional image: the vortex lines, connecting the centres of the elementary cell faces, 
are drawn as solid lines |62| . 



magnetic field can be added in the Hamiltonian |l| by substituting — 0j with 0^ — (pj + Aij , where Aij = A ■ dl 
and the magnetic field B is given by ,B = rotA. Naturally, the lines tend to have the same winding direction. In the 
case of neutral particles it must be noticed that the description in terms of a magnetic field is given in the rotating 
frame, so that the vortices are actually rotating in the inertial reference frame. 

We also investigated the case of the coupling ( 16 1, that is derived from the NS superfluid gap. By fitting the results 
of [IS] with the function [A + Bz) e~'-^^ , we see that the gap is spatially asymmetric with respect to the centre of 
the crust. The parameter B controls the ratio of the slopes of the gap at the inner and outer crust boundaries. From 
the fit we get that this ratio is « 5. This means that the inner vortex wall (bottom side in the figure) is much thinner 
than the outer wall and it has a sharp interface with the vortex-free region as shown in Fig. [7] 

Regarding the application of the Hamiltonian ([I]) to rotating NS, we observe that using the inhomogeneity (16) in 
presence of a strong magnetic field gives results qualitatively similar to those of Fig. [6] 

The correlation function (cos (0o — 0r)) between the site at the origin and the site at distance r is plotted in Fig. 
|8j one clearly sees that the correlations goes to zero in a length scale <^ determined from the condition (3J{£,) ^ 0.5 . 
In the same figure we plot the average number of vortices at distance r from the centre, which goes from to a finite 
value at the same length scale. 

To compare our results with LDA, we calculated the average number of vortices for a homogeneous system at 
various temperatures and compared it to the local density of vortices as a function of the coupling strength in Fig. [9] 
The LDA agrees well with the inhomogeneous results, so it is expected that local observables are accurately described 
by LDA. 

Fig. |8] clearly illustrates that the number of vortices n„(r) is large in the region of small couplings (where the 



Molecular Physics 



8 




FIG. 2: A typical configuration above the critical temperature in the homogeneous limit (/3J — 0.4). Here the vortices are not 
separated and the inter-vortex distance is smaller than their radius. 



correlation with the central region is vanishing): therefore we introduce an effective number, He// given by 

n,ff{r)^l3J{r)n4r). (17) 

The Fig. 10 plots the effective number of vortices ( 171 for the same parameters of Fig. [s] illustrating the anticipated 
fact that the effective number of vortices is maximum in the region with intermediate couplings. 



CONCLUSIONS 



In this paper we studied the distribution of vortex lines in inhomogeneous three-dimensional lattices modelled by 
an XY model with spatially varying coupling constants. The motivation of our work was two-fold and related to the 
superfluid properties of rotating trapped ultracold atoms in optical lattices and rotating neutron stars. Indeed the 
low-energy properties of ultracold atoms in deep optical lattices are in general described by Hubbard-like models: in 
the superfluid phase, either due to temperature and/or to the effect of a rotation, vortex lines will develop in the 
system. Since an external potential is usually present in order to trap these systems, then an arising issue is the 
determination of how the trapping potential affects the distribution vortex lines. 

The study the effect of inhomogeneities on the vortex lines in superfluids in a lattice is also relevant in the context of 
neutron stars. Theoretical calculations and indirect observational evidence support the idea that crustal free neutrons 
form a superfluid. Since the inhomogeneous neutron superfluid in the inner crust has a gap exhibiting a strong spatial 
dependence and it is experiencing the periodic potential created by the lattice of very neutron rich nuclei, then an 
interesting problem is the study of the effect of the inhomogeneity on the vortex lines in the inner crust. Also relevant 
would possibly be the individuation of a simplified model capturing the main properties of vortex lines in neutron 
stars. 
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We have presented results for the spatial distribution of the number of vortex lines, showing how the vortex lines 
are expelled from the region with higher couplings. The reason for this behaviour is that in the homogeneous limit the 
critical temperature is proportional to J and for temperatures smaller than a few small vortex loops are thermally 
generated. In presence of inhomogeneity, there is a space dependent value of the temperature value at which the vortex 
loops start to blow out; thus, with the couplings varying across the lattice, in the regions with small coupling constants 
vortex loops accumulate, while in contrast vortices are expelled from the regions with large couplings. Defining an 
effective number of vortices to take into account that where the coupling is small the number of superfluid particles 
is also small, one has that the effective number of vortices will reach a maximum value at intermediate values of the 
couplings. Results for the correlation functions and in presence of a strong external magnetic field (or strong rotation 
for neutral superfluids) are also presented, 

As a future work, it would be interesting to determine the vortex-vortex interaction in presence of inhomogeneities 
and/or dislocations: from this point it would be useful to write the inhomogeneous Hamiltonian in terms of dual 
variables and write the vortex-vortex interactions starting from the (inhomogeneous) lattice Green's functions. The 
corresponding results would then allow for the determination of the effective size of the region with an appreciable 
effective number of vortices. Equally interesting would also be the systematic study of the effects of an (eventually 
frustrated) external magnetic field. 

In perspective, we think that a promising line of work is the extension of the results obtained in this paper using a 
simplified XY toy model to quantitatively study vortex behaviour in rotating neutron star crusts. Attempts in this 
direction have been made using a phenomenological Gross-Pitaevskii equation 63 . The planned extension consists 
in writing an effective Hamiltonian for the superfluid neutrons put in rotation and feeling an effective inhomogeneous 
periodic potential: the superfluid dynamics will be then described by a simplified (ATy-like) model having parameters 
determined from the underlying microscopic fermionic Hamiltonian. The resulting model could be used to make 
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FIG. 4: The vortex configuration for the spherically symmetric inhomogeneity (14 1 with R = 3 and /3Jo = 1. Vortices are 
shown as solid lines where the coupling J is maximum and vanish proportionally to the coupling. Where /3J (r) < 0.1 the lines 
become practically not visible. There is a central region of size R where the vortex lines do not penetrate, and a vortex shell 
around this core. 



predictions for observable quantities and comparisons with available data, with the aim to improve the current state 
of the art both in terms of modelling and of input parameters. Using numerical results from microphysical simulations 
as inputs, for the pairing gap, the effective lattice potential and the inter-particle interaction, one could be then able 
to study the behaviour of vortices in a realistic way and hopefully bridge the gap between the astrophysical vast 
dataset and the many-body theoretical picture describing the crust. 
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FIG. 8: Numerical results in the spherical geometry (Eq. (14 1 with 13 Jo = 1, -R = 4, 6). The size of the lattice is L = 21 and 



the distances are measured in units of the lattice spacing. The correlation between the site at the origin r — and the one at 
distance r, C(r) = {cos(0(O) — (pir))) is plotted for 7? = 4 (x) and R = 6 (A). The average number of vortices n is also plotted 
(□ for R — 4, Q for i? = 6) to show how vortices screen the correlation. 
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FIG. 9: The number of vortices as a function of the local temperature for the same model as Fig. [s] (□ for R — A and Q for 
R = &) and for several homogeneous systems (A for LDA). 
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FIG. 10: The effective number of vortices at radius r for the same spherical geometry and parameters as in Fig. [s] (i? = 6). 



